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Abstract. We develop a new thermodynamic approach to stochastic graph-rewriting. 
The ingredients are a finite set of reversible graph-rewriting rules Q (called generating 
rules), a finite set of connected graphs V (called energy patterns), and an energy cost 
function e which associates real values to each of these energy patterns. The idea is that Q 
defines the qualitative dynamics, by showing which transformations are possible, while V 
and e allow one to attach an energy to the reachable graphs and, thereby, describe their 
long-term probability distribution tt. Given Q and V, we construct a finite set of rules Qv 
which (i) has the same qualitative transition system as Q-, and (ii) when equipped with 
rates according to e, defines a continuous-time Markov chain of which n is the unique fixed 
point. The construction relies on the use of site graphs and a technique of ‘growth policy’ 
for quantitative rule refinement which is of independent interest. 

This division of labour between the qualitative and long-term quantitative aspects of 
the dynamics leads to intuitive and concise descriptions for realistic models (see It also 
guarantees thermodynamical consistency {aka detailed balance), otherwise known to be 
undecidable, which is important for some applications. Finally, it leads to parsimonious 
parameterizations of models, again an important point in some applications. 


1. Introduction 

Along with Petri nets, communicating finite state machines, and process algebras, models 
of concurrent systems based on graphs and graph transformations (GTS) have long been 
investigated as means to describe, verify and synthesize distributed systems [13] . Beyond 
their visual aspect, which is often useful in modelling situations, there is a lot to like about 
GTSs: there are category-theoretic frameworks to express them and encapsulate their syntax; 
and the existence of a strong meta-theory m is a reassurance that methodologies developed 
in specific cases can be ‘ported’ to other variants. 

2012 ACM CCS: [Software and its engineering]: Software notations and tools—Eormal language 
definitions—Syntax / Semantics; [Theory of Computation]: Logic; Formal languages and automata theory— 
Formalisms—Rewrite systems. 
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Graph-rewriting rules are convenient for writing compact models and modifying them [8] , 
and lend themselves naturally to probabilistic extensions [I8l|20]. However, for all their 
flexibility, even rules can only do so much. We ask in this paper “what if we did not have 
to write the rules?”. This is where we take a page from the book of classical statistical 
mechanics. In such models, which often involve graph-like structures, as in the Ising model, 
the dynamics is not described upfront. Instead, the system of interest is equipped with an 
‘energy lansdcape’ which specifies its long run behaviour, be it deterministic as in classical 
mechanics, or probabilistic in statistical physics. The dynamics just follows from the energy 
data. In the eye of a computer scientist, this use of energy looks like a latent syntax. (This 
is especially true in the application of these ideas to molecular dynamics.) 

The broad aim of this paper is to make this syntax explicit by introducing energy patterns 
and costs from which the total energy of a state of the system can be computed; and to dehne 
a procedure whereby, indeed, the dynamics described as probabilistic graph-rewriting rules 
can be derived from these energy data. Descriptively, this takes us to an entirely new level of 
conciseness (as in the example shown in Q. It also guarantees thermodynamical consistency 
by construction, a property known otherwise to be undecidable m- This property provides 
a predictable equilibrium state for the system which can be used as a guide when writing 
models by giving an intuitive understanding of favoured states and the (possibly overlapping) 
subgraphs within them. But perhaps the nicest byproduct of this approach is the fact 
that the methodology leads to parsimonious parameterizations. The parameter space which 
usually scales as the number of rules (which in turn has at best a logarithmic impact on the 
cost of a simulation event 0), will now scale as the number of energy patterns provided in 
the specification. 

The particular GTSs we consider form a reversible subset of the Kappa site-graph 
stochastic rewriting language. Kappa is used for the simulation and analysis of combinatorial 
dynamical systems as typically found in cellular signalling networks |22[ |2^ and has been 
predicted to “become one of the mainstream modelling tools of systems biology within the 
coming decade” [1]. Similar graph formalisms where nodes have a controlled valence/degree 
have been considered e.g. the BNG language miiH], Kissinger and Dixon’s quantum proof 
language m, and Kirchner et al. chemical calculi [3] . Site-graph rewriting has found recently 
a ‘home’ both in the single-pushout GTS tradition [9] and the double-pushout one ISIEI- 
This makes one hopeful that the thermodynamic methodology we propose can cross over 
to other fields where quantitative GTSs can be used, e.g. in the modelling of adaptive 
networks [16]. While our scalable energy-based parameterization is particularly important in 
biological applications where parameters often need to be inferred, one can imagine it to be 
useful in other modelling situations with uncertainty. 

Outline: We start by introducing a running example that will help us illustrate the concepts 
presented throughout the paper. Then we proceed with the dehnition and relevant properties 
of the specific GTS we use, namely a simple reversible fragment of Kappa. Next, we introduce 
growth policies (adapted from Ref. [25]), a tool which allows one to replace a rule with 
an orthogonal set of refined rules while preserving the quantitative semantics. We use this 
tool with a specific policy which refines a rule into finitely many rules, each of which has 
a definite energy balance with respect to a given set of energy patterns. This leads to our 
main theorem which guarantees that the stochastic dynamics of the obtained refined rule 
set converges to an equilibrium distribution parametrized by the cost of each energy pattern. 
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Throughout, the presentation is set in category-theoretic terms and mostly self-contained. A 
substantial example concludes the paper. 


1.1. Running example: Assembling triangles. The following example will be used to 
show the intrincacies of our rule generation mechanism. We consider graphs with three types 
of nodes where each can only bind nodes of a different type, and at most one thereof. As 
generators, we consider simple binding and unbinding rules subject to this constraint (see 
details in { 2.4): 


0 • 



• o 


m-o 

o ® 


0-# 


These rules, given an unbounded supply of the three types of nodes, can create chains 
(of any length) and closed cycles (of length some multiple of three) as shown in Fig. Given 
the simplicity of the graphs that can be formed, we can describe them using a linear textual 
notation where numbers are used to represent the three types of nodes, then chains are 
simply written as words {e.g. 2312) while for cycles we indicate half-edges at both ends, e.g. 
•123- is the triangle and •123123- is the hexagon. We will use that shorthand notation in the 
sequel. 

Our goal here is to be able to favour the formation at equilibrium of certain structures 
like the triangle by assigning a significantly negative energy to them (the convention is that 
a lower energy means a higher probability). 

2. Site graph rewriting 


2.1. Site graphs and homomorphisms. A site graph G consists of finite sets of agents 
(nodes) and sites (connection slots), Ac and So, a partial function ac '■ Sq Ac-, and a 
symmetric edge relation Eg on Sq- The pair Sq, Eg form an undirected graph; sites not in 
the domain of Eg are said to be free. The role of the additional map cjg is to assign sites to 
agents; sites not in the domain of og are said to be dangling, and will be used to represent 
half-edges (see below). Usually one also endows agents and/or sites with states (as we do in 
the example treated in Q; our main construction in ^carries over readily to these. 

One says G is realizable iff (i) sites have at most one incident edge; (ii) no dangling site 
is free; and, (iii) edges have at most one dangling site. Note that point (i) implies that no 
site has an edge to itself. 


A homomorphism h \ G ^ G' oi site graphs is a pair of 
functions, hs ■ Sg —)• Sg' and h_A_ '■ Ag —^ Ag', such that (i) 
whenever hj\^[aG{s)) is defined, so is crG'{hs{s)) and they are ac 

equal; and (ii) if s Eg s' then hs{s) Eg' hs{s'). 

A homomorphism /i : G —?■ G' is an embedding iff (i) hj[ Ag ^ 

and hs are injective; and (ii) if s is free in G, so is hs{s) in 
Gh If /i : G —)• G' is an embedding and G' is realizable then G is also realizable. 


Sg' 

Ag' 
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Site graphs and homomorphisms form a category SG with the natural ‘tiered’ composi¬ 
tion; embeddings form a subcategory; if in addition, we restrict objects to be realizable, we 
get the subcategory rSGe of realizable site graphs and embeddings. 

Returning to our chains-and-cycles example, we can use sites as a means to enforce on 
nodes our earlier constraint where a node can connect to at most one node of each other 
type. Here we present a possible encoding of a simple chain of 3 nodes: 



Ag = {x, y,z} , ac = {sa ^ a : s € {l,r} ,a € Ac} , 

•Sg = Xx, ly, Vy, Iz, V z"\ , £g — {(Txi ly): {}yi X x), {Xy, Iz), {Iz^ ^J/)} • 

Note that this site graph is realizable and the lx and sites are free. To complete the 
encoding of our example we need types for nodes which we now introduce. 


2.2. The category of site graphs over C. A homomorphism /i : G —>• (7 is a contact 
map over C iff (i) G is realizable, (ii) ac is total and (hi) whenever hs{si) = hs{s 2 ) and 
c’'g(si) = o'g('S 2 ); then si = S 2 - The third condition of local injectivity means that every 
agent of G has at most one copy of each site of its corresponding agent in C; C is called the 
contact graph. 

Hereafter, we work in the (comma) category rSGec whose 
objects are contact maps over C, and arrows are embeddings 
such that the associated triangle commutes in SG. We write 
Tc{h, h') for the set of such embeddings between /i, h' contact 
maps over C; we also write |_| for the domain functor from 
rSGec to rSGe which forgets types. In particular, if h : G — )• C 
is a contact map, we write \h\ for its domain G. 

Note that in rSGe, an embedding h : G ^ G' may map a dangling site s of G to any 
site s' of G' (provided hs is injective). In particular, s' does not have itself to be dangling. 
This means that dangling sites can be used as an “any site” wild card when matching G. In 
the typed case, i.e. in rSGec, the contact map h : G ^ C tells us which agent in G the 
dangling site belongs to because uc is total, and this must be respected by h. Dangling sites 
can now be used as “binding types” to express the property of being bound to the site s of 
an agent of a given type. 

The contact graph G is fixed and plays the role of a type: it specifies the kinds of agents 
that exist, the sites that they may possess, and which of the l^cKl^d + l)/2 possible edge 
types are actually valid. It also gives canonical names to the types of agents and their sites. 

In our running example we use the triangle G = •123- as contact graph: 




•^c = {1)2)3} , Sc = {h,ri,l2,r2,l3,r3} , 

(Tc = {sa*-^ a : s € {l,r} ,a € Ac} , 

£c = {{ri,l2), {l2,ri), {r2,l3), {l3,r2), {r3,h), {h,r3)} 
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This contact graph constrains agents to bind only agents of a different type. The local 
injectivity condition on contact maps restricts the number of sites available for each node to 
two. Our simple chain of three nodes from the previous section ( ^2.1[ ) can thus be typed 
using a contact map h : G ^ C with /i _4 = {x i-A 1, ?/ 1 —)• 2, 2 ; i—>■ 3} and hs = {s^; i-a • 

s £ {l,r} ,x £ Ag}- 


We now extend the shorthand notation introduced in 11.1 to cover the case when hs is 
not locally surjective. Specifically, when a node does not mention its two allowed sites, we 
write ‘?’ for the missing site. A missing site can be introduced by an embedding and can be 
free or bound. Sites do not explicitly show in this notation but they are implicitly positioned 
at the left {1) and right (r) of each agent. Binding types are denoted by an exponent. Thus, 
for instance, in ^12? we have two agents u, v of type 1, 2, respectively, where u is bound to a 
dangling site of type r^, on its lu site and bound to u’s ly site on its Vy site. Similarly, v is 
bound to u on its ly site and does not have any r site. Using this extended notation we can 
write down our three rules in text: ?1 + 2? ^ ?12?, ?2 + 3? ^ ?23?, and ?3 + 1? ^ ?31?. 


Recall that the goal of this example is to control the type of cycles which one finds at 
equilibrium. To this effect, we introduce a set of energy patterns V consisting of one contact 
map for each edge type (ie. ?12?, ?23?, and ?31?), and one for the triangle •123-. As we will 
see, a careful choice of energy costs for each pattern will indeed lead to a system producing 
almost only triangles in the long run (see Fig. in Appendix [A]) . 


2.3. Multi-sums in rSGec. The category SG has all pull-backs, constructed from those 
in Set; it is easy to see that they restrict to rSGe^. SG also has all push-outs and all sums, 
but these do not in general restrict to rSGec. (Just as push-outs and sums in Set do not 
restrict to the subcategory of injective functions.) 



Figure 1: An example of a mixture of chains and 3n-gons which our three rules can reach 
starting with 12 agents of each type. 
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However, rSGec has multi-sums meaning that, for 
all pairs of site graphs of type C, hi : Gi ^ C and 
/i 2 : G 2 —> C, there exists a family of co-spans 9\ : hi ^ 

Si /i 2 : @ 2 ) such that any co-span 71 : /ii —)• h /i 2 : 

72 factors through exactly one of the family and does 
so uniquely. In outline, given a co-span 71 : hi —)• h <— 

/i 2 : 72 , we first take its pull-back (which is guaranteed 
to remain within rSGec) and then the push-out of that: the fact that the original co-span 
exists implies that this push-out also remains within rSGec. The multi-sum is then a choice 
of push-out for each (isomorphism class of) pull-back. 

The idea is that the pairs hi, enumerate all minimal ways in which one can glue hi 
and /i 2 , that is to say all the minimal gluings of Gi and G 2 that respect G. There are finitely 
many, all of which factor through the standard sum in the larger slice category SGc- 

The notion of multi-sum dates back to Ref. |12) : we will call them minimal gluings 
in rSGec according to their intuition in this concrete context, and use them in §3.3| to 
construct balanced rules. 

To illustrate this idea, let us consider, in the context of our example, the minimal gluings 
of Di = 7123? and D 2 = 72317. Computing them is a matter of computing all pull-backs. 
One such pull-back is Di •(— 7237 —)■ D 2 , which gives us the minimal gluing 712317. On the 
other hand, the span Di ^ 737 —)• D 2 is not a pull-back since for all cospans that can close 
the square, the span factors through Di ^ 7237 —)■ D 2 . This is a consequence of 737 not 
being a maximal overlap of Di and D 2 '. whenever 3 is contained in a possible overlap of Di 
and D 2 , 2 will also be there. 

All minimal gluings of Di and D 2 are displayed in the following diagram with their 
respective pullbacks. The left square has the empty graph as pullback and 71237 -|- 72317 as 
pushout, the central square has 7237 and 712317, and the right square has 7237 -|- 717 and 
•123-. Each square uses arrows of different colour and style. 


7 02 , 

hi -^ Si ^- h2 
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2.4. Rules. A rule r over C is a pair of contact maps : 
L^C,rn:R^C which differ only in their edge structures, 
i.e. Al = Ar, Sl = Sr, ai = (Tr, = rR,A and vr^s = 
rR,s- 

For example, the unbinding rule gi 2 =?12? —)• ?1 + 2? 
may be represented as: 



Ar 

= Ar 

= {n,u} 

Sr 

= Sr 

- ^v} 

ctr 

= CFr 

= {ru>-^ u,ly v} 

rL,A 

= rR^A 

= {u 1—>■ 1, u 1—>■ 2} 

rL,s 

= rR^s 

= ri,ly i-A I2} 


with only the edge structures Sr = {{r^, Iv), {Iv, ^u)} and Sr = 0 differing. 

A contact map /i : G —>■ C is a mixture iff ac is total (no dangling edge) and hs is 
locally surjective, i.e. for all a G Ac-, hs{aQ^{a)) = a^^{hA{ci)). In words, a mixture is a 
fully-specified site graph with respect to the type C. 

An embedding : ur ^ h induces a rewrite of 
the mixture h by modifying the edge structure of 
the image of from that of to that of rR. The 
result of rewriting is a new mixture h*, where \ h*\ has 
the same agents and sites as \h\, and an embedding 
Ip* : rR^ h*. 


rR 


tr 


h* 


( 2 . 1 ) 


This type of rewriting can be formalized using double push-out rewriting [H], but with 
the simple rules considered here, there is no need. The inverse of r, defined as r* := [ rr , rr ) 
is also a valid rule; by rewriting h* with r* via ip*, we recover h and fi. 


Given a finite set of rules G over C, we define a labelled transition system Cg on mixtures 
over C: a transition from a mixture /i is a rewriting step determined and labelled by an 
event {r,ip), as in diagram (2.1), with r in ^ and ip in Tc{r, h). We suppose hereafter that 
G is closed under rule inversion, i.e. G = G*. Hence, every (r,i/’)-transition has an inverse 
{r*,ip*), and Cg is symmetric. 


2.5. CTMC semantics. It is not difficult to see that for any rule r, \ Tc{rR, h)\ < 
where d(r) is the number of connected components in vr . Hence, Cg has finite out-degree, 
bounded by \G\ • for some d. Also, as agents are preserved by rules, the (strongly) 

connected components of Cg are finite. Given a rate map k from G to M-|-, we can therefore 
equip Cg with the structure of a time-homogeneous irreducible continuous-time Markov 
chain (GTMG), simply by assigning rate k{r) to an event of the form (r, ip). We write Cg 
for the GTMG thus obtained. 

A finite time-homogeneous GTMG Ai has detailed balance for a probability distribution 
TT on Al’s state space iff, for all states x and y, 7r(x) • q{x, y) = 7r(y) • q{y, x) where q{x, y) 
is A4’s transition rate from x to y. This implies, assuming Ai is irreducible, that vr is the 
unique fixed point of the action of Ai to which the probabilistic state of Ai converges, 
regardless of the initial state. 

In our case, the probability 7r(x) will be proportional to where E{x) is the energy 

of the system as defined by a set of patterns and their associated costs. 
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2.6. Extensions and rule refinement. We say an embedding (j) 
(/>' : s —)• s" if there is some embedding 6 ■. s' ^ s" such that 
9(j) = (j)' and write ip < ^ for this. We refer to a prefix of an epi 
(/> : s —)• as an extension of s. In the category of extensions of 

s, a morphism between objects <p : s ^ s' and cp' : s ^ s" is an 
embedding 6 : s' ^ s" such that the triangle commutes. If 0 is an 
iso we write (p=s p' ■ 


s —)> s' is a prefix of 



Epis of rSGec are characterized as follows |2S]- suppose s : G ^ C and s' : G' ^ C 
are contact maps then (p : s ^ s' oi rSGec is an epi iff every connected component of G' 
contains at least one agent in the image of Rule application preserves epis and in fact 
also prefixes of epis; 


Lemma 2.1. Let r = {rL,rR} be a rule and p : ri ^ x be an embedding with ri, xr, and x 
contact maps in rSGec. The embedding p* : xr ^ x* that results from applying r to p is a 
prefix of an epi iff p is. 


Proof. This amounts to proving that ip* > p* is an epi p > p is. For this it is enough 
to consider the case where the rule adds or deletes exactly one edge. Rules that modify more 
than one edge at a time can be decomposed as sequences of deletions and insertions of edges. 
Given that each deletion and insertion preserves the property, the sequence will preserve it as 
well. The case of adding an edge is easy, as the image of p* has fewer connected components 
to “touch”. The case where r deletes an edge can introduce new connected components, 
however in this case both ends u, v of the deleted edge must be in r^, so whether the 
deletion disconnects or not the codomain of p, the components of p*{u) and p*{v) will have 
a pre-image, namely u and v. □ 

It follows that the category of extensions of xr and xr are isomorphic. Indeed, because 
we have chosen to restrict to reversible rules, the full categories under xr and xr are evidently 
isomorphic, and by the lemma above this correspondence preserves being prefix-of-an-epi; if 
p is an extension of xr, we will write p* for the corresponding extension of xr. 

A family of epis pi : s ^ ti uniquely decomposes s, or is a refinement of s, if, for all 
mixtures h and embeddings p : s ^ h, there exists a unique i and p' such that p = p'pi. 
This is the basic requirement for a reasonable notion of rule refinement: it guarantees that 
the LHS s of a given rule splits into a non-overlapping and exhaustive collection of more 
specific cases ti. 

In the next section, we will be constructing specific such decompositions in order to 
produce families of subrules that are compatible with energy patterns, i.e. each subrule 
should produce and consume the same number of energy patterns in each application, 
regardless of the particular embedding p G Tc{rR, h) that the subrule is applied to. 

To find such unique decompositions we first recall the growth policy method [25], which 
works by detailing which agents and sites should be added to s. Specifically, a growth policy 
r for s is a family of functions indexed by all extensions p : s ^ t, where maps 
u G A\t\ to 8- subset Pcpiu) of {tj\^{u)), i.e. each agent in |f| is allocated a subset of the 
set of sites belonging to the agent t_ 4 ^{u) it is mapped to in G. An agent in |t| may cover 
some, or all, of these sites or even completely extraneous sites: if the former, i.e. for all 
u in A\t\, (u)) C rfp{u), we say that p is immature; if for all us, the inclusion is an 

equality and p is also an epi, we say that p is mature; otherwise p is said to be overgrown. 
The functions must satisfy, for all extensions p and p' > p, the faithfulness property. 
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= ^<j>' ° i’A with V’ such that ijjil) = 4>'', so a site requested by cj) must be requested by any 
further extension. Additionally, this property forces F to eagerly ask for all sites that will 
be eventually requested at any given agent in the codomain of (j). If (p is not overgrown then 
no (j)' < (j) is overgrown either. Also, note that the union of two growth policies is itself a 
growth policy. 

Given an s and a growth policy F for s, we define F{s) by choosing one representative 
per =s-isomorphism class of the set of all extensions of s which are mature according to F. 

The following theorem guarantees that factorizations through F{s) are unique when 
they exist, but not that they necessarily do exist. In the next section, we will construct a 
specific growth policy of interest for which this property of exhaustivity of the decomposition 
can be proved by hand. As such, it fulfils our desired criteria of providing an exhaustive 
collection of mutually exclusive sub-cases. 

Theorem 2.2. Let s and x he contact maps and F a growth policy for s. If an embedding 
Ip : s ^ X in Tc{s,x) can be decomposed in two ways as •jipi and ^2(p2 with (pi : s ^ F in 
F{s) and 7* : — >■ x then pi = p2 and 71 =72. 


Proof. Suppose that 'jipi = F 2 P 2 , where 
pi and p 2 are mature extensions of s accord¬ 
ing to F and pi 7 ^ p 2 . As shown in dia¬ 
gram ( 2 . 2 ), we have an inner square formed 
by the pull-back tti, 712 , and the minimal 
gluing 01 , 02 of 71 and 72 . Also 61 and O 2 
are epis, as every connected component of 
m has a pre-image in ti or t 2 and so also in 
s, since the piS are epis, and so also in the 
other of t 2 and ti. 



( 2 . 2 ) 


If 01 , 02 are not both isomorphisms then there must be a pair u, z, consisting of a node 
in m with pre-images ui, U2 in ti, t2 and a site 2 ; of u, such that z has no pre-image in 
exactly one of 0i, 02 - Say it is 62- Since pi is not overgrown, 2 ; G T'(^j(ui) and, by faithfulness, 
2 ; G F(j,{{ui,U 2 )), where {ui,U 2 ) is the pull-back pre-image of ui and U 2 . So again, by 
faithfulness, 2 ; G ^<^ 2 (^ 2 ) which contradicts our original assumption. So 0i and O 2 are isos. It 
follows that pi = p2 as there is only one representative per =s-isomorphism class in F{s). 
Finally, 71 = 72 because pi is an epi. □ 

Given a rule r and an extension p : ^ t of ri, we write for the refined rule 

associated to p] that is, is the pair (t, F) with t* the codomain of the extension p* 
corresponding to p. Given F a growth policy for ri, we write F{r) for the family of rules 
obtained by refining r according to F ; that is F{r) is the family of rules r^, for p ranging in 

r{rL)- 

An example of growth policy is the ground policy which assigns all possible sites to all 
agents. In this case, F{s) is simply the set, possibly infinite, of all epis of s into mixtures, 
considered up to =s', and F{r), the ground refinement of r, contains all refinements of r 
along those epis, which therefore directly manipulate mixtures. It is easy to see that the 
ground refinement for our running example is infinite, since each of the three rules can 
trigger the extension of a chain of any length. 
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3. Rule generation 


We fix a finite set G of generator rules; and a finite set V of connected contact maps in 
rSGec; these are our energy patterns. 

A canonical set of generator rules can be derived from the contact graph C by constructing 
a pair of binding/unbinding rules (and a full set of state changes if we were to use site graphs 
with states as we do in applications) for each edge in £c- This set is maximally general in 
that each generator rule asks for the least possible context to trigger a binding or unbinding 
event. This is the set which we have chosen to use in our running example. In the next 
example which we present in Q we will consider a strict subset of this canonical generating 
set. But, in general, there is no prescription regarding which rules one decides to incorporate 
to G- 


The goal is now to refine G into a new rule set Gv where each refined rule is V-balanced.! 
which means that, however applied, it consumes or produces a fixed amount of each c in 
V. The construction proceeds in two steps. Firstly, we characterize balanced refinements; 


secondly, we define a specific growth policy with balanced mature extensions. Using Th. 2.2 
we show that these mature extensions obtain a proper finite refinement of G. 

Note that ground extensions of g are trivially balanced but, in general, the ground 
refinement is impractically large or even infinite; ours will always be finite. 


3.1. Consumption and production of patterns. Consider c in "P and a rule r. For an 
r-event to consume an instance 7 of c in a mixture h, 75, and "05 must have images which 
intersect on at least one site which is modified by r (by adding an edge if it was free, or 
by deleting or changing an edge it was part of). This is the 
case iff the associated minimal gluing ( 7 ',!/^') (obtained by 
restricting the co-span to the union of its images in h) has 
the same property. Likewise, for an r-event to produce an 
instance of c, the associated minimal gluing between c and 
rji must have a modified intersection. We call such minimal 
gluings relevant] they are the ones which underlie events that 
can affect the instances of c. 

This notion can be illustrated by looking at the minimal gluings of the left-hand side of 
the unbinding rule ?12? —)• ?1 -|- 2? and energy pattern 71231? (supposing for a moment that 
we had 71231? in V): 
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Of these, only the central square is a relevant minimal gluing since the image of the edge 
consumed by the rule is contained in the image of the energy pattern in the minimal gluing. 



(3.1) 


3.2. 'P-balanced extensions. Given g vn. Q with LHS ql and 4> : gi ^ t an extension of 
gL, we say that (j) is V-left-balanced iff, for all relevant minimal gluings 'j : c ^ m t : 9 
with c € V, 9 IS an isomorphism. This means that the image of c under 7 is contained in t. 
Symmetrically, one says that (j) is V-right-balanced iff (p* 
is a P-left-balanced extension of r*. An extension cp is 
V-balanced iff it is P-left- and "P-right-balanced. We say 
that a balanced extension cp is prime iff it is minimally so, 
i.e. any prefix of (p that is P-balanced is isomorphic to (p 
(as an extension). Prime extensions are epis since erasing 
an ‘untouched’ connected component in the codomain 
preserves balance. 

If (/> is a P-balanced extension of g, the refined rule g^ has a balance vector in Z^, 
written Acp, where, for each c G P, A(p{c) is the number of copies of c produced by 
any 5 ^-event, which is also the difference between the number of embeddings of c in the 
RHS and the LHS of 5 ^. In other words, for any mixture h, balance of cp guarantees 
A(p{c) = I Tc{c,h*)\ -\Tc{c,h)\ = \ Tc(c, fl'</,,R)| - | Pc(c,5?i,T)l- Conversely, if c in P violates 
the condition of diagram (3.1), well-chosen applications of g^f, will result in different and 
context-dependent A<p{c). Thus, the notion of balanced extension characterizes the property 
that we want. 


To illustrate these definitions consider the unbinding rule gi 2 =?12? —?1 -|- 2? with 
the trivial extension (p = 1^^ (the identity arrow of the left-hand side of rule <712 )• Define 
Pi = {?12?, ?23?, ?31?} and P 2 = {T23-}. On the one hand, <p is Pi-balanced. Of all possible 
gluings of patterns in Pi with (p^s codomain, only the trivial one of ?12? with itself is relevant 
to < 712 . On the other hand, cp is not P 2 -balanced, because there is a gluing between •123- and 
?12? mapping the 12 edge to itself in the triangle. This gluing is relevant since applying 7712 
breaks the triangle; and clearly 9 : ?12? —?■ •123- is not an iso. Failure to be balanced reflects 
in the fact that different applications of 7712 incur different values of A(p{-12‘i-)-. namely 1 or 
0 depending on whether the broken edge belongs to a triangle. 


The ideas introduced now are discussed in a detailed example in (3.5 


3.3. Absorb or Avoid. We have just seen how the property of being P-balanced charac¬ 
terizes, among all extensions cp ol gi, those with an unambiguous energy balance. The next 
step is to define a growth policy with a set of mature extensions which are balanced and 
form a unique decomposition of 52 ,. If we return to diagram (3.1) in the case where 9 is not 
an isomorphism, the natural idea is to request the addition to the codomain of cp of these 
sites which are in m and belong to nodes in the image of 0^. In this way, we force cp to grow 
further and make a decision about whether it wants to absorb the image of c in m, or would 
rather avoid this by growing in a way that is incompatible with the 7 , 9 gluing. 


Some care is needed to ensure faithfulness, i.e. P^ = F^/^cp'^, since relevant minimal 
gluings on cp can disappear along a further extension cp' so that a site that was requested at 
cp may no longer be so after at cp'cp. To address this, we add site requests from all relevant 
minimal gluings in the past of an extension. 
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Given 51 in ^ we define a growth policy 
r for gi. Suppose cj) : ql ^ t is an exten¬ 
sion of ql- We set to request a site s in 
cr^^(tyl(u)) at agent u in A\t\ either (i) 
u = 4>Aiuo) and s = 4>s{so) for some uq in 
So in 5 |g^|; or (ii) (p factorizes as 4>2pi, 
where (pi ■ 9l ^ ti, and there is a relevant minimal gluing 7 : c — )• m' •(— : 9, with c in V, 

and some ui in and a site s'^ in ( 0 yi(ui)) such that u = (p2,A{'^i) and s = m^(s'^). 

We refer to the extension 02 : ti —)• t as a rewind of 0; we say that the request of s at tt 
originates from ui. The first clause simply ensures that all sites already covered in gi are 
asked for; the second one adds in sites which appear by gluing at some point between gi 
and t, and implements the absorb-or-avoid constraint explained beforehand. 

Symmetrically, we dehne a growth policy F* for gji by applying the same definition to 
the reverse generator g*. Since extensions of gL and gn are isomorphic, we can, with a slight 
abuse of notation, define F^ := F U F*. 

Theorem 3.1. The above F^ is indeed a growth policy for gi^; the induced refined family of 
rules F^{g) is exhaustive, non-empty, V-balanced, and finite. 



Proof. We take the same notations as in diagram (3.2). 


Growth policy: Clearly, T'^(ui) C F^{u) as every request for a site in ti will propagate 
to t by definition. To prove the other direction, we need to verify that the requests generated 
by rewinds do not depend on the choice of factorization. So, wlog, consider gluings on left 
extensions of g, and let an alternative factorization of 0 through t2 be given which gives rise 
to a site request in u originating from some U 2 in t 2 : 



t 3 u 




Consider p the pull-back of the two rewinds (ie the lower co-span); by construction it 
contains a pre-image for ui and U 2 , say {ui,U 2 ). The relevant minimal gluing of c and t 2 
that makes the site request restricts to another (minimal) gluing of c and p. This new gluing 
is still obviously relevant (as it contains the same overlap with the original gi)] as such, the 
same site request is made by the pre-image (ui, U2) agent in p which then propagates to ui 
in ti as required. 


Surjectivity: Let p : gL ^ x he an embedding of gi into a mixture. We can restrict the 
co-domain of 0 to be the connected closure y of the image of 0 in x, resulting in an epi 
Py '■ 9l ^ V- Let us further restrict y by removing: 1) all unneeded agents, i.e. those that 
have no sites requested by the growth policy, 2 ) all sites not requested by the growth policy 
except those bound to sites that are requested and, finally, 3) all connected components that 
no longer have a pre-image in gi', call the result z. The non-requested sites left in 2 : act as 
binding types (see 12.2). We thus obtain an epi pz : gL ^ z which is mature with respect to 
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rv 

since, by construction, it has all the sites requested by (j)y and so, by faithfulness, all 
those requested by (pz- 

Non-empty: Clause (i) guarantees that we request at least the sites in g which implies 
that g is not overgrown. 

Balanced: If (/> G B'^{g) is not balanced then there must be some relevant minimal gluing 
inducing a further site request, hence (p cannot be mature. 

Finite: A request for a site a at some node in an extension cp : gi ^ t, oi (p* : gji ^ t, 
originates from a relevant minimal gluing of some c in "P with a prefix (pi of (p. Because this 
gluing is relevant, it must be that a is at a distance from the image of gL in the codomain of 
(pi which is at most (5(c), the diameter of c (else c would not intersect the image of gi). The 
same bound holds in the codomain of (p, as distances can only contract by further extension. 
Therefore any site requested in t has a distance to the image <p{gL) which is bounded by 
maxcgp (5(c). If (p is not overgrown, this sets a bound on the diameter of t. Hence there are 
finitely many mature extensions. □ 

Therefore, given Q and V, we obtain a finite P-balanced rule set which rehnes Q 
exhaustively, by setting Qp := (jg^gF^(g) (disjoint sum). To every refinement g^, corresponds 
an inverse refinement g'^*; hence, G-p = Qj, is closed under inversion like Q. 

3.4. Other options considered. As said earlier, there are other ways to generate a 
balanced rule set. Beyond the ground extensions, one idea is to use primes, i.e. minimal 
balanced extensions. Clearly, prime extensions factorize all balanced ones -including the 
ones we generate by using the absorb-or-avoid growth policy presented above. The problem 
is that this factorization is not unique: prime extensions do not dehne in general a valid 
refinement as distinct sub-cases may overlap. 

Another idea is to obtain balanced rules by gluing V directly onto a generator rule g, 
in all possible maximal combinations, rather than on all extensions of g as in the above 
growth policy. This always reveals enough of the context in which g is applied so as to get 
P-balanced refinements because, by definition, no additional form can be glued further. The 
problem with this approach is the opposite to that with primes: it does not in general build 
enough rules; so the refinement is valid but not exhaustive. 

3.5. Triangles in detail. We now discuss the concepts introduced in this section using the 
example of the unbinding rule gi 2 = ?12? —)■ ?1 -f 2?. Of the four patterns in V, we need 
only consider the triangle T = •123- as the others all trivially glue or trivially don’t. 

We apply our growth policy to the extension 
(pi7 : ?12? —7- 12? whereby we reveal a new free site 
in 1. The relevant minimal gluing T —)■ T 712? 
disappears in this extension, as the new site is 
bound in T. However, by rewinding back to 7127, 
our growth policy requires us to reveal additionally 
the second site of 2; as such, (pi? is not mature 
despite being balanced (indeed minimally so). 

Intuitively, the growth policy forces us to reveal additional sites of 712? because an 
instance of gi 2 may or may not consume a triangle and so has an ambiguous energy balance. 


712? 
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The definition of the growth policy makes it explicit that we have to reveal both of the 
hidden sites of ?12? and so there are four sub-cases to consider: 12, ^12, 12^ and ^12^ (using 
the notation introduced earlier for binding types). The first three cases all avoid the triangle, 
thereby having definite energy balance, and are moreover mature with respect to the growth 
policy; however, the case of ^12^ requires further analysis as it covers both the case of a 
triangle and of any chain of length at least 4. 

Consider then the case of the extension to 
73123? where the two binding types of ^12^ have 
been embodied in two distinct 3 agents. Any 
rewind of interest has to go back far enough to 
glue T; there are two maximal such extensions: 

: 712? ^ 7312? and 07,37 : 712? ^ 71237. 

The former requests the hidden site of the left- 
hand 3 agent of 731237; the latter (diagram not shown) requests the hidden site of the 
right-hand 3 agent. As such, 73123? is not mature with respect to the growth policy and we 
must expand it into its four sub-cases: 3123,^3123,3123^,^3123^. 

Our final refinement, according to the growth policy, consists in the following extensions 
12, ^12, 12^ •123-, 3123, ^3123, 3123\ ^3123^ 

If were to glue only on 7127, not on all extensions thereof, we would only generate the rule 
T —)■ 231, i.e. the case where the bond happens to belong to a triangle. Clearly, this is a valid 
refinement but is far from covering all cases. Note also that the extension 0 i ,2 : 712? —?■ 12 
is not a minimal balanced extension as it factors through either 01 ,? : 712? —?■ 12 ? or 
07,2 : 712? —)• 712 which are both minimal. This illustrates the problem with prime extensions 
discussed earlier: the candidate refinement {127, 712} is not valid precisely because it leads 
to an ambiguous decomposition of the case of 0 i ,2 and will therefore generate rates that are 
different for instances of gi 2 with the same balance vector A0. 

On the other hand, the subrules 3123, ^3123, 3123^, and ^3123^ all share the same 
balance vector. In fact, they can be summarized by a single prime extension 073,3? ■ ?12? —>• 
73123? which, together with 0^ : 712? —?• T, mutually exclusively and exhaustively refines 
03x23 : 712? —7- ^12^. Our growth policy is too local to see this optimisation, and has to refine 
0 ? 3 , 3 ? into the above four explicit subcases. 

That said, the output of the growth policy could be either optimized by hand or indeed 
subjected to the procedure of rule compression which would automatically perform this 
optimization for us. A working Kappa model where the generators of this example have 
been fully expanded according to the growth policy and then manually compressed can be 
found in Appendix [A| 



73123? 


3.6. Rates. To equip Gp with rates, we suppose given a "P-indexed real-valued vector of 
energy costs e, and a rate map k : Qp ^ M+ such that, for all g^ in Qp: 

^og k{gl^) -log k{g^) = e ■ Acf) (3.3) 

with A0 in Z^, the balance vector of the refined rule g^ with respect to P, a well-defined 
quantity by Th. |3.1[ 

We write P(x) for the P-indexed vector which maps c to | Tc{c,x)\, and define the 
energy E{x) of x as e • P(x), i.e. the sum over all c G P of the product of the number of 
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instances of c in x and the energy cost of c. We also write Cg{x) for the hnite (strongly) 
connected component of x in Cg, and define a probability distribution (in Boltzmann format) 
TT^ on Cg{x) by: 


TTx{y) ■= 

We can now prove our main theorem. 


^-e-Viy) 




(3.4) 


Theorem 3.2. Let Q, V, Q-p, k, and tTx he defined as above; then (i) Cg^, and Cg are 
isomorphic as symmetric LTSs; and (ii) for any mixture x, the irreducible time-homogeneous 
continuous-time Markov chain Cfg^ has detailed balance for, and converges to, on Cg^{x). 


Proof. Both Cg and Cg.p offer transitions from a mixture x: the former are labelled by 
pairs with g in Q, in Tc{gL,x); the latter by pairs [g,f>,^) with g^ the refinement of 

g along a mature extension cp : gp ^ t, and 7 in Tc{t, x). Steps in the latter can be mapped 
to steps in the former by transforming labels as follows: (5^,7) e->■ [g,'y(p). As Gp refines G 
exhaustively (Th. 3 . 1 ), this correspondence is a bijection, which establishes the first claim. 

(Pedantically, there is a full and faithful functor between the two corresponding free 
categories which is the identity on objects —incidentally, this bijection is readily seen to 
respect the symmetries on labels.) 


Since we have multiple rules in Cg^,, each of which can be applied in several ways, there 
can be more than one transition from x to the same y —each uniquely described by a 
(5^,7) label. Each such {grfi,'y) has an inverse, (fl'^*,7*), where: g* is the rule inverse to g; 
(p* corresponds to (p in the isomorphism between the categories of extensions of x and y, 
with (p_A = 7* is the embedding corresponding to 7, also with 7 a = 7^- One can 

easily verify that (p* is an epi, and that (p* is also mature. Hence (5^*, 7*) determines a valid 
transition in Cg^, which is inverse to (5,/,, 7), and we have a bijective correspondence between 
transitions from x to y and those from y to x. 


Consider a pair e, e* of such corresponding events due to g^ and 


P*; because 


e IS a 


transition from x to y, and (p is P-balanced (Th. 3 . 1 1, we have V{y) = 'P(x) + Acp, and hence 
e • Acp = e ■ {V{y) — V{x))\ so, by ( | 3 . 3 [ ), the rates of e, e* are such that: 

A:(e*) e“'^'^ 0 ) = k{e) 


and by summing this equation over all pairs, we obtain detailed balance for the probability 
local to the component Cg.^{x) = Cg.p{y), defined above as -Kx = 'Xy, since: 

q{y,x) = q{x,y) 

The convergence statement follows by standard results on finite CTMCs. □ 

Note that the subset of the state space which is reachable from x in Cg, namely 
Cg{x), is hnite; hence, the partition function Z{x) := &Cg{x) ® ^ which hgures in the 

denominator of tTx is also hnite. In the presence of rules which increase the number of agents, 
the components Cg{x) can be inhnite and Z{x) may diverge. For (mass action stochastic) 
Petri nets, convergence is guaranteed if detailed balance holds, but it is not true in general 

for Kappa mm- 

Another point worth making is that the result holds symbolically —regardless of the 
energy costs e. So e can be seen as a set of parameters, an ideal support for machine learning 
techniques if one were contemplating htting a model to data. 
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3.7. A linear kinetic model. We now ask: what of the actual rates of Among all 
possible choices which accord with (3.3), we delineate a tractable subset whose size grows 
quadratically in \V\. This is a useful log-linear heuristic, which is common in machine 
learning, but has no particular claim to validity. 

Suppose we have, for every generating rule in a constant Cg G M, and a matrix Ag 
of dimension \V\ x \V\. Subject to the constraints that Cg* = Cg, and Ag* + Ag = I, we can 
dehne a log-affine rate map which satisfies (3.3) by: 

^og{k{g^)) := Cg - Ag(e) • Aip (3.5) 

The kinetic model expressed in (3.5) requires of the order of iPp x |^| parameters. In practice, 
one needs even fewer parameters, as only those energy patterns that are relevant to a given 
g, i.e. have non-zero balance for at least one rule in r^{g), need to be considered when 
building Ag. Typically, for larger models, this will be a far smaller number than \V\. This 
relative parsimony is compounded by the fact that the number of independent parameters 
will be often lower, because the Acf) family often has low rank. It is to be compared with the 
total number of choices which is far greater as it is of the order of the number of refinements, 
that is to say l■^^(5)l• 

If we set Cg* = Cg = 0, Ag* = 0, Ag = I, we get: k{g^) = k{g^i,) = 1. As 

€■ A(p is the difference of energy between the target and source in any application this 
choice amounts to being exponentially reluctant to climb up the energy gradient. This is a 
continuous-time version of the celebrated Metropolis algorithm |24j . 

Another particular case, completely symmetric and which can be a reasonable choice to 
begin with, is obtained for Ag* = Ag = 1/2: 

kg{4>) 

kg* {(/A = 

with Cg = As A(jA = —Acj), this is indeed a symmetric definition. 

Finally, it is fun to draw a comparison between the ascription given in (3.5) and the 
Arrhenius rate law. This law posits a dependency of the rate constant /c of a reaction of the 
form log/c = c — Ea/kT, where c is a constant (defining the basic time scale of the reaction), 
Ea is the so-called activation energy of the reaction and T is the temperature. In our case, 
we are not concerned with the effect of T on the (logarithm of the) rate but with the effect 
of consuming and producing various energy patterns in V at the locus of the instance of 
the generator rule g. In this view of things, (3.5) posits that the ‘activation energy’ of (p 
depends linearly on the cost of the various patterns and the balance of p. 


= Cge -^-^‘^/2 


3.8. Energy functions do not need to be linear. We now return to a key assumption 
made in the preceding section and consider a more general situation where the energy 
function E is no longer asked to be linear. For reasons to become clear shortly, we still 
assume the much weaker property that E can be factored as v oV{-) for some finite set of 
patterns V. 

rSGec—^ mSet(P)—M (3.6) 

Note that if we see multi-sets as equipped with the usual point-wise partial order, 'P(_) is 
evidently functorial. 









THERMODYNAMIC GRAPH-REWRITING 


17 


As an example, consider a divalent agent with contact graph X(a^,b^) —where the 
shared superscript symbolizes an edge between sites a and b — with which one can form only 
chains and cycles, and a single generator g which can create/delete the unique edge type. 
Write C 3 for a cycle of length 3 (a triangle), and t2 for an open chain of length 2. We define 
the quadratic energy function E{x) = \ J’c(c 3 ,x)p, i.e. v{n) = n^. Applying g forward to a 
chain t2 in a site graph of the form x = t 2 + x' will create a new copy of C3, and give the 
following energy balance: 

AE = {\Tcic3,x)\ + l)‘^ - {\Tcic3,x)\f = 2| rc(c 3 , x)| + 1 

Therefore, detailed balance forces the log-ratio Kg of the backward/forward rates assigned 
to an edge creation to depend on x. This is unlike the case of linear potentials where this 
ratio is independent of x. However, this extension g^p of g (where it is applied to a chain 
t2) is balanced with respect to V. This means, as we have seen, that the stoichiometric 
"P-vector Acj) associated to g^ —where each component A(f>{c) is defined as the difference of 
I ^c(c, 2/)| — I ^c(c, x)| for a ^'-transition from x to y — is the same for all x, y. In the example 
Acj) has only one component, which is constant A4>(c3) = 1 because the binding of the two 
free extremes of an open chain t2 can only produce one triangle, regardless of the context in 
which the refined rule g^ is applied. 

In general, one can visualize the situation as follows. 

9 

X - ^y 

V 


V 


V 


w— 


+AE 


and detailed balance amounts to asking for Kg = v{V{x) + A4>{x)) — v{V{x)). If v happens 
to be linear then this is the usual condition Kg = v{A(f){x)). If v is not linear, the constraint 
does not seem very helpful as a priori one has to know x to compute the right hand side. 
But by the assumption (3.6) opening the paragraph, Acf factors through P, hence: 


Kg = v{j:>(^x)+ag{V{x)))-v{V{x)) =: ?AgOp(x) 

where the second equation dehning ifg uniquely as a real-valued function from P-multisets. 
With this rewriting, it is plain that the constraint depends not on full knowledge of x, but 
only on P(x). Equivalently, we see that Kg factors through P(-) just like E. In the example, 
kg = 2| Tc{c 3 , x)| -I- 1, and ipgin) = 2n + l. 

This is good enough to define rates for a balanced g^. For example, by analogy with the 
earlier linear kinetic model, we can choose log-rates (seen as real-valued multi-set functions) 
as follows: 


kg — Olg (^•'^) 

with Og, j3g real-valued functions on P-multisets such that Og* = Ug and fig* fig = 1. This 
assignment solves the above constraint as, clearly, fjg* + ftg = 0. 
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From the simulation point of view, this added generality requires two things: (i) that rates 
can be made to depend explicitly on observables; (ii) that the internal state of the simulation 
be extended to incorporate V{x). Both possibilities are already generically available in the 
current version of the main Kappa simulator KaSim [2]. A slight modification of the engine 
(not implemented) could obtain direct updates to 'P(x) as, by assumption, applying leads 
to a constant +A(/) update; and the same holds for propagating these updates to the rates 
of the rules which depend on them, e.g. as in (3.7). Thus, the complexity properties of the 
simulation algorithm are preserved m- 


4. Allosteric ring 

We now put our methodology to use on a realistic example of a bacterial flagellar engine. In 
this section, we will use the traditional syntax of Kappa to denote site graphs: subscripts 
for states and shared superscripts for edges between sites, e.g. A(xq), Unlike the 

mathematical definitions of ^ agent and site types are indicated as explicit labels. Again, 
we use KaSim (the standard Kappa engine) for the simulation shown below. 


4.1. The model. The engine can rotate clockwise or anti-clockwise at high angular velocities, 
and this will decide whether the bacterium tumbles or swims forward. One can build a 
simple model of the switch between the two modes [2]. The engine is seen as a ring of n 
identical components, P, with two possible conformations, 0 and 1. (In reality, each of the 
n = 34 component protomers is itself a tiny complex made of different subcomponents, but 
the model ignores this.) A ring homogeneously in state 0 (1) rotates (anti-) clockwise and 
induces tumbling (straight motion). Importantly, neighbouring Ps on the ring prefer to have 
matching conformations. States of the ring with many mismatches thus incur high penalties. 
In the absence of any Y molecule binding a P, its favoured conformation is 0; conversely, 
in the presence of a Y, P favours 1. {Y stands for a small diffusible protein named CheY.) 
To bind, Y has to be activated by an external signal. Hence the switch can be triggered by 
a sudden activation of Y which then binds the ring and induces a change of regime. The 
sharper the transition between the two regimes the better. 

As each of the Ps can be in four states, the ring has on the order of I0^° non-isomorphic 
configurations which precludes any reaction-based {e.g. Petri nets) approach to the dynamics 
where each global state is considered as one chemical species. At this stage, we could apply 
the rule-based approach, or, better, we can obtain the rules indirectly by applying the 
methodology of ^ This is what we now do informally. 

First, we define our contact graph with two agent types: P{x, y, /o,i, s) with domains x, 
y to form the ring, s to bind its signal Y, and / a placeholder for P's conformation; Y{su,p) 
with two internal states to represent activity. 


Second, we capture the informal statements in 
the discussion above by defining the energy patterns 
and associated costs. Note that the various motifs 
overlap. Following ^ we associate to each ring 
configuration x the occurrence vector V{x) and 
total energy e • V{x). For example, a ring of size 
n uniformly in state 0 and with no bound Y s has 


Motif 


Cost 

P{fi,x^ 

P{f^) 

),P{y\fj) 

ePP 

P{h,s^ 

),y{s^) 
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total energy n{eQQ + e^). This, in turn, defines the equilibrium distribution of the ring, 
namely x has probability proportional to exp(—e • V{x)). (The convention is that the lower 
the energy, the likelier the state.) 

In order to complete our energy landscape, we need 
to pick energy costs which reward or penalize local con¬ 
figurations as discussed above: the role of (4.1) is to 


pp pp 

mo ) ^11 

^PY 


< 

< 

> 


ef 


pp 


align the internal states of neighbours on the ring —an 
Ising penalty term for mismatching neighbours which will 
“spread conformation”; (4.2) makes 0 the favoured state, 

while (4.3), which kicks in only in the presence of T, makes 1 the favoured state. 
The next step is to create the dynamics. The naive rule b for PY binding: 

b:=P{s),Y{sp)PYP{s^),Y{sl) 


01 ( 4 . 1 ) 

(4.2) 

(4.3) 


has a AE which is ambiguous as it will be either Cq 


PY 


or e 


PY 


depending on its instances; 


hence, we have no hope of assi gnin g rates to this rule that satisfy detailed balance —unless 
= ef^, which contradicts (4.3). To get a definite balance, one needs to refine this rule: 

bo := P{fo, s),Yisp) ^ P{fo, s^),Y{sl) 
bi := P(/i, s),Y{sp) ^ P{h,s^),Y{sl) 

Now each rule specifies enough of the context in which it applies to have a definite energy 
balance Following the same intuition of revealing (just) enough context, we obtain a 
balanced rule set for conformational changes: 

fij ■= Pix^, , s), Pix"^, fj) o P{fi,y^),P{x^Ji,y‘^,s),P{x^,fj) 

f'ij ■= P{fi,y^),P{x^,fo,y'^,s-),P{x‘^, fj) ^ P(/i,y^),P(x^/l,y^s-),P(x^/J) 

The first (second) group of rules represents the changes in the absence (presence) of a T 
bound to the middle P undergoing a change of conformation. (The fact that P’s site s is 
bound is indicated by the underscore exponent.) 

These /-rules have respective balance: 

pp , pp _ pp _ pp ,^P_^P 
fPP _L .PP _ .PP _ .PP <.P_.P< .PY _ .PY 

As we have ten reversible rules, and only eight energy patterns, there must be linear 
dependencies between the various balances. Indeed, in this case, it is easy to see that the 
family of vector balances has rank six. Thermodynamic consistency induces relationships 
between rates; a well-established fact in the case of reaction networks {e.g. see Ref. Cl). 

With the rules in place, the final step is to choose rates which satisfy detailed balance. 
This guarantees that the obtained rule set converges to the equilibrium specified by the 
choice of the energy cost vector. Convergence will happen whatever e is, i.e. symbolically. 
If, in addition, e follows (4.1-4.3), one can see in Fig. [^that the ring (i) undergoes sharp 
transitions when active Y is stepped up and down again; and (ii) has at all times very few 
mismatches. 
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Figure 2: The simulation steps up the amount of active Y aX t = 100, and down again at 
t = 200; this sends the entire ring into an homogeneously 1 conformation, and back 
to 0. The number of mismatches (lowest curve) stays low, even during transitions. 


4.2. How to generate the rules. Our set of balanced rules was based on two generators, 
h for binding, / for flipping: 

b:=P{s),Y{sp)^P{s^),Y{sl) 
f := P(/o) o P(/i) 

Note that there is a design choice here. In effect, we are saying that we are not interested in 
forming/breaking the bonds between the Ps in the ring. If we wanted to incorporate also 
the ring assembly in the model, we would have to add P{x),P{y) -H- P{x^), P{y^) among 
our generator set Q. This would generate many more rehned rules, as we will see. Recall 
that our patterns fall in three subgroups: P{fi,x^),P{y^, fj)] P{fi)] and P{fi, s^),Y{s^). 

Consider the extensions of b: clearly only the last pattern can glue relevantly on it; the 
corresponding (unique) site request is for P to reveal / and its internal state. This gives the 
first two rules bo, bi. 

Consider now the more interesting extensions of /: the second pattern type glues 
relevantly but does not generate any site request; the third one asks P to reveal its site s, 
resulting in two possible extensions {s- means that s is bound): 

P{fo,s)^P{h,s) 

P(/o,s-)oP(/i,s-) 

These extensions are not mature yet, as one can glue relevantly patterns of the first type on 
both sides of P, inducing a further request for revealing P’s sites x and y. 

If we are in the component of an initial state where Ps are arranged in a ring, then we 
know that the neighbours on both sides exist and are Ps; this gives the hnal rehnement of 
the above into the rules /y, /b described earlier. If, on the other hand, we do not know that. 
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Figure 3: Snapshots of the ring configuration are taken at time 500, 1500, and 2500. Solid 
(green) circles indicate conformation 1, hollow ones conformation 0; a dot in the 
centre indicates a bound (hence active) Y. At times 500, 2500, no Y is bound 
(because they are all inactive) and the ring is globally in state 0, up to tiny 
fluctnations; at time 1500, it is globally in state 1 as a consequence of the binding 
of ys. 


we also have to add several rnles where one or both of x, y are free, corresponding to open 
P-chains. This demonstrates the sensitivity of the obtained rule set to the initial choice of 
generators. 


Hence, the rules we generated by hand are indeed the ones we would generate using our 
general refinement strategy of ^ 

We can visualize the obtained simulations by extracting snapshots before, after and 
during the injection of active Ys, as in Fig. Again we see few mismatches in both regime 
because of the Ising interaction expressed by the energy costs. The choice of rates made 


m 


Ref. [2] for the /-generator is the symmetric version of (3.5). 


5. Conclusion 

We have presented a new ‘energy-oriented’ methodology for the development of site graph 
rewriting models based on a set V of energy patterns; these patterns use a graphical syntax 
which allows us to specify the energy landscape. Rewrite rules are implicitly defined by V 
and generator rules Q. The resulting rule set Qp is guaranteed to be thermodynamically 
correct and to converge to the probability distribution described by the energy landscape, 
given suitable rates. The construction is entirely parametric in the energy costs e, and 
modular in Q. This means that in a modelling context, one can sweep over various values for 
e without having to rebuild the model, and compositionally add new rule components to G. 
Both features are clearly useful. We expect our construction to provide a broad and uniform 
langnage to describe and analyse models of interacting biomolecules and similar systems of 
a qnantitative fine-grained and distributed nature. 

There are no specific conditions bearing on this construction other than that energy 
patterns should be local. It would be interesting to investigate whether suitable constraints 
on patterns and generator rules can obtain optimized generated rule sets. Another interesting 
extension would be to deal with non-local forms of energies expressing long-range interactions, 
where the metric is read off the graph itself. In practice, there will be many more rules 
generated, and beyond the descriptive aspects, simulations will need new ideas to be feasible. 
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A ray of hope comes from the log-affine kinetic model presented in §3, as rules can be 
partitioned by energy balances for faster selection. 

Finally, as said in the introduction, there is a growing body of literature which turns 
a theoretical eye to site graph rewriting [ElIlslIISlE], and it is tempting to ask whether 
our derivation can be replayed in more abstract settings; in particular, it would be very 
interesting to investigate its integration with the abstract framework for rule-based modelling 
developed in [23]. 
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Appendix A. Triangles model 


A working model definition for our first example (introduced in 11.1) is spelled out in full 
here. The syntax used is for KaSim version 4. As in the example of §4, we use the symmetric 


linear kinetic model of 13.7 


A.l. Agents, Parameters and Initial State. We start by declaring agents together with 
their sites. Then we declare the energy costs for T = •123- (‘t’), ?12? (‘ab’), ?23? (‘be’), ?31? 
(‘ca’) as parameters. Subsequently the initial state is defined as containing 1000 agents of 
each type. 

# Agent signatures 
7oagent: A(l,r) 

7oagent: B(l,r) 

7agent: C(l,r) 

# Energy costs 
7var: ’t’ -10 
7var: ’ab’ 0 
7var: ’be’ 0 
7var: ’ca’ 0 

# Initial state 

7init: 1000 (A(), B(), C()) 

A.2. Observables. Next we declare the graphs we would like to keep track of, i.e. those 
graphs we would like to generate trajectories for plotting. Here we are interested in the 
fraction of agents that are used to assemble triangles. 

7obs: ’T’ |A(1!1, r!2), B(l!2, r!3), C(l!3, r!l)| 


A.3. Rules. Finally, we write our set of rules. These have been grouped according to the 
qualitative generator rule they refine. Note that these rules have been manually compressed 


as explained in 13.5 


# A(r), B(l) -> A(r!l), B(l!l) refines into: 

A(l,r), B(l,r) -> A(l,r!l), B(l!l,r) @ [exp] (-1/2 * ’ab’) 
A(l!r.C,r), B(l,r) -> A(l!r.C,r!1), B(l!l,r) <3 [exp] (-1/2 * ’ab’) 
A(l,r), B(l,r!l.C) -> A(l,r!l), B(l!l,r!l.C) @ [exp] (-1/2 * ’ab’) 
A(l!l,r ), B(1 ,r!3), C(l!3,r!l) -> \ 

A(l!l,r!2), B(l!2,r!3), C(l!3,r!l) @ [exp] (-1/2 * (’ab’ + ’t’)) 
C(r!l), A(l!l,r ), B(1 ,r!3), C(l!3) -> \ 

C(r!l), A(l!l,r!2), B(l!2,r!3), C(l!3) @ [exp] (-1/2 * ’ab’) 


# A(r!l), B(l!l) -> A(r), B(l) refines into: 

A(l,r!l), B(l!l,r) -> A(l,r), B(l,r) @ [exp] -(-1/2 * ’ab’) 
A(l!r.C,r!l), B(l!l,r) -> A(l!r.C,r), B(l,r) @ [exp] -(-1/2 * ’ab’) 
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A(l,r!l), B(l!l,r!l.C) -> A(l,r), B(l,r!l.C) @ [exp] -(-1/2 * ’ab’) 
A(l!l,r!2), B(l!2,r!3), C(l!3,r!l) -> \ 

A(l!l,r ), B(1 ,r!3), C(l!3,r!l) 0 [exp] -(-1/2 * (>ab' + ’t’)) 

C(r!l), A(l!l,r!2), B(l!2,r!3), C(l!3) -> \ 

C(r!l), A(l!l,r ), B(1 ,r!3), C(l!3) @ [exp] -(-1/2 * ’ab’) 

# B(r), C(l) -> B(r!l), C(l!l) refines into: 

B(l,r), C(l,r) -> B(l,r!l), C(l!l,r) 0 [exp] (-1/2 * ’be’) 

B(l!r.A,r), C(l,r) -> B(l!r.A,r!1), C(l!l,r) 0 [exp] (-1/2 * ’be’) 

B(l,r), C(l,r!l.A) -> B(l,r!l), C(l!l,r!l.A) 0 [exp] (-1/2 * ’be’) 

B(l!l,r ), C(1 ,r!3), A(l!3,r!l) -> \ 

B(l!l,r!2), C(l!2,r!3), A(l!3,r!l) @ [exp] (-1/2 * (’be’ + ’t’)) 

A(r!l), B(l!l,r ), C(1 ,r!3), A(l!3) -> \ 

A(r!l), B(l!l,r!2), C(l!2,r!3), A(l!3) ® [exp] (-1/2 * ’be’) 

# B(r!l), C(l!l) -> B(r), C(l) refines into: 

B(l,r!l), C(l!l,r) -> B(l,r), C(l,r) 0 [exp] -(-1/2 * ’be’) 

B(l!r.A,r!l), C(l!l,r) -> B(l!r.A,r), C(l,r) 0 [exp] -(-1/2 * ’be’) 

B(l,r!l), C(l!l,r!l.A) -> B(l,r), C(l,r!l.A) 0 [exp] -(-1/2 * ’be’) 
B(l!l,r!2), C(l!2,r!3), A(l!3,r!l) -> \ 

B(l!l,r ), C(1 ,r!3), A(l!3,r!l) @ [exp] -(-1/2 * (’be’ + ’t’)) 

A(r!l), B(l!l,r!2), C(l!2,r!3), A(l!3) -> \ 

A(r!l), B(l!l,r ), C(1 ,r!3), A(l!3) ® [exp] -(-1/2 * ’be’) 

# C(r), A(l) -> C(r!l), A(l!l) refines into: 

C(l,r), A(l,r) -> C(l,r!l), A(l!l,r) 0 [exp] (-1/2 * ’ea’) 

C(l!r.B,r), A(l,r) -> C(l!r.B,r!1), A(l!l,r) 0 [exp] (-1/2 * ’ca’) 

C(l,r), A(l,r!l.B) -> C(l,r!l), A(l!l,r!l.B) 0 [exp] (-1/2 * ’ca’) 

C(l!l,r ), A(1 ,r!3), B(l!3,r!l) -> \ 

C(l!l,r!2), A(l!2,r!3), B(l!3,r!l) @ [exp] (-1/2 * (’ea’ + ’t’)) 

B(r!l), C(l!l,r ), A(1 ,r!3), B(l!3) -> \ 

B(r!l), C(l!l,r!2), A(l!2,r!3), B(l!3) @ [exp] (-1/2 * ’ca’) 

# C(r!l), A(l!l) -> C(r), A(l) refines into: 

C(l,r!l), A(l!l,r) -> C(l,r), A(l,r) 0 [exp] -(-1/2 * ’ca’) 

C(l!r.B,r!l), A(l!l,r) -> C(l!r.B,r), A(l,r) 0 [exp] -(-1/2 * ’ca’) 

C(l,r!l), A(l!l,r!l.B) -> C(l,r), A(l,r!l.B) 0 [exp] -(-1/2 * ’ca’) 
C(l!l,r!2), A(l!2,r!3), B(l!3,r!l) -> \ 

C(l!l,r ), A(1 ,r!3), B(l!3,r!l) @ [exp] -(-1/2 * (’ea’ + ’t’)) 

B(r!l), C(l!l,r!2), A(l!2,r!3), B(l!3) -> \ 

B(r!l), C(l!l,r ), A(1 ,r!3), B(l!3) <S [exp] -(-1/2 * ’ea’) 

A.4. Results. We run this model using KaSim (version 4) to obtain trajectories for the 
number of triangles during the simulation. From these numbers we can estimate what is 
the expected number of triangles at equilibrium. In Fig. [^the results of 4 different runs 
with different energy costs for a unit triangle are displayed. When e(T) = —10, almost all 
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Time (arbitrary units) 


Figure 4: Trajectories for T = •123- when e(T) varies (this value is set by changing the value 
of parameter ‘t’ in the KaSim file). 


agents are used to build triangles. Instead, when e(T) = —2.5 only less than 5% is used. 
Interestingly, in both cases the set of states that minimize the energy function is the same, 
namely those states that maximize the amount of triangles. So then why is it that in the 
latter case there are so few triangles? The reason is entropic: although the probability 
of being in a state with few triangles is small, there are many such states and together 
they outweigh the probability of being in the few states were the energy is minimized. By 
further decreasing the energy of those few states we compensate for this mass effect, until at 
e(T) = —10, order wins, and the effect is not noticeable anymore. 
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